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ABSTRACT 

HELAC is a FORTRAN based package that is able to compute efficiently helicity 
amplitudes for arbitrary scattering processes within the standard electroweak 
theory. The algorithm exploits the virtues of the Dyson-Schwinger equations 
as compared to the traditional Feynman graph approach. All electroweak ver- 
tices are included in both the unitary and Feynman gauges, and computations 
including all mass effects are available. A version performing multi-precision 
computations with arbitrary - user defined - accuracy is also included, allowing 
access to any phase space point for arbitrary high energies. 
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PROGRAM SUMMARY 



Title of the program: 
HELAC. 

Catalogue number: 
Program obtainable from: 

Dr. Costas G. Papadopoulos, Institute of Nuclear Physics, NRCPS 'Democritos', 15310 

Athens, Greece, e-mail Costas .PapadopoulosOcern. ch. 

Also at ftp : / /alice . nuclear . demokritos . gr/dist/pub/ papadopo/helac . 



Licensing provisions: none 

Computer for which the program is designed and others on which it has been tested: 
ALPHA, HP, IBM workstations. 

Operating system under which the program has been tested: 
Unix. 

Programming language: 
FORTRAN 77 and FORTRAN 90 

Keywords: 

Automatic evaluation of helicity amplitudes, Dyson- Schwinger equations, recursive algo- 
rithms. 

Nature of physical problem: 

A substantial part of particle phenomenology nowadays is based upon our ability to study 
efficiently processes involving a relatively large number of particles. This requires efficient 
algorithms for matrix element calculation and phase space generation and integration. 
As far as the matrix element calculation is concerned, the traditional approach utilizes 
standard Feynman graph representation of the scattering amplitude, resulting to a compu- 
tational cost that grows asymptotically as n!, where n is the number of particles involved 
in the process. 

Method of solution: 

As an alternative recursive algorithms based on Schwinger-Dyson equations lead asymp- 
totically to a much lower growth of the computational cost, namely a", where a ~ 3. It is 
the aim of this paper to present such an algorithm as well as the corresponding FORTRAN 
code which allows the computation of any tree-order electroweak amplitude. 
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1 Introduction 



A substantial part of current particle phenomenology is based upon our ability to study 
efficiently processes involving a relatively large number of particles. This requires efficient 
algorithms for matrix element calculation and phase space generation and integration. 
As far as the matrix element calculation is concerned, the traditional approach utilizes 
standard Feynman graph representation of the scattering amplitude, resulting to a compu- 
tational cost that grows asymptotically as n\, where n is the number of particles involved 
in the process. As an alternative recursive algorithms based on Schwinger-Dyson equa- 
tions lead asymptotically to a much lower growth of the computational cost, namely a", 
where a ~ 3. It is the aim of this paper to present such an algorithm as well as the 
corresponding FORTRAN code which allows the computation of any tree-order electroweak 
amplitude. 



As advertised, the algorithm is based upon the well-known Dyson-Schwinger equations. 
These equations are equivalent to the field equations derived from the classical lagrangian 
and express recursively the n-point Green's functions in terms of the 1—, 2—, . . . , (n — 1)- 
point functions. 

In order to illustrate this point let us consider a simple example of a QED-like theory 
possessing minimal interactions of a spinor field to a gauge boson. Let pi,p2, ■ ■ ■ ,Pn 
represent the external momenta involved in the scattering process under consideration, 
taken to be incoming. For a vector field we define 



a four vector, which describes any sub-amplitude from which a boson V with momentum 
P can be constructed. The momentum P is given as a sum of momenta of external 
particles, namely P — J2ieiPi with / C {1, . . . , n}. Accordingly, we define by 



a four-dimensional spinor, which describes any sub-amplitude from which a fermion with 
momentum P can be constructed and by 



2 The algorithm 






(2) 




(3) 
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a four-dimensional antispinor, which describes any sub-amphtude from which an an- 
tifermion with momentum P can be constructed. 

The Dyson-Schwinger equations take the following simple form. For a boson with 
momentum P'*, 



+ 



i=l P=Pi+P2 





where the propagator 11 is given by 



p2 _ ^2 



and e(Pi, P2) is a sign factor, that takes the value ±1 and whose definition will be described 
below. For a fermion with momentum P 





m Y.^p=pMPi)+ E it9)rp{P2MP,)e{p,,P2) 

P=Pl+P2 



1=1 



where V is given by 



f-m 



and for an antifermion with momentum P, 



+ 





where 



(5) 



^(P) = Y.5p=,MP^)+ E {^gmPlW{P2)Pe{P,,P2) (6) 

i=l P=Pi+P2 



The scattering amplitude can be calculated by any of the following relations, 

{bQ{Pi)b^{pi) where i corresponds to a photon 
il)Q{Pi)il){pi) where i corresponds to an incoming fermion line (7) 
'ijj{pi)'ijjo{Pi) where i corresponds to an outgoing fermion line 

where 

Pi = Y.Pj^ 

so that Pi + Pi = 0. The functions with subscript are given by the previous expressions, 
except for the propagator term. This is because the outgoing momentum Pj being on shell 
the propagator factor is removed by the amputation procedure. The initial conditions are 
given by 



b''ip^) = e^fe),A = ±l,0 

r uxiPi) ifp°>0 
\ vx{-Pi) ifp°<0 

r uxiPi) ifp°>o 

1 vxi-Pi) ifp° <0 



where the explicit form of e'^,ux,vx,ux, vx are given in the Appendix. 

In order to actually solve the recursive equations it is convenient to use a binary 
representation of the momenta involved For a process involving n external particles 
with momenta pf , i = 1 . . . ,n we assign to the momentum P'^ defined as 

p^=j:pf 

where / C {l,...,n}, a binary vector m = (mi, . . . , m„), where its components take the 
values or 1 in such a way that 



P^ = Y.mi pf. 

1=1 

Moreover this binary vector can be uniquely represented by the integer 

n 
i=l 

with < m < 2"^^. Especially, the external momenta Pi, i = 1 . . . ,n are represented by 
m = 2*~^, i = 1, . . . ,n. All momenta P can now be replaced by the corresponding integers 
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A very convenient ordering of integers in binary representation relies on the notion of 
level /, defined simply as 



i=l 



As it is easily seen all external momenta are of level 1, whereas the total amplitude 
corresponds to the unique level n integer 2"~^. This ordering dictates the natural path of 
the computation; starting with level- 1 sub-amplitudes, we compute the level-2 ones using 
the Dyson- Schwinger equations and so on up to the final expression Eq.(^. As far as the 
sign factor is concerned 

e(Pi,P2) e(mi,m2) 

we define 

e(mi,m2) = (-l)^('"^''"^) (9) 

with 

ma) = Y.^^i\Yl 1 (10) 




where hated components are set to if the corresponding external particle is a boson. It 
is not difficult to see that this sign factor takes properly into account the anti-symmetry 
of the amplitude with respect to fermionic particles. 

As an illuminating example we present here the computation of the amplitude for the 
process e^e" The calculation starts with the computation of the level-2 sub- 

amplitudes, which are all sub-amplitudes produced from any two of the external particles. 
These are 

6^(12) = (z(7)n^2^^(4)7.^(8) 
^(18) = {ig)^{2)P{lQ)P,s 
Vi(20) = iig)^iAmiQ)V2o 
^-(24) = (^(7)^24^16)^(8) 

Then we compute the level-3 sub-amplitudes: 

^(14) = (z^7)V^(2)Kl2)Pi4 

6^(28) = (z(7)n^8^(^(20)7.V^(8)+V^(4)7,V^(24)) 
and then the level-4 sub-amplitude, which is the final level in this case, 

^o(30) = (tg) {i^{2)m) + ^(14)^16) + ^(18)^(12)) (11) 
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Then the amphtude is simply given by 



A = V^o(30)V(l) 

Note that we have chosen the particle number 1 as our ending point, so we have computed 
all sub-amplitudes where the momentum pi does not appear: this excludes all odd integers 
between 1 and 2""^. 

There are two special issues in the computation which go beyond the recursive equa- 
tions presented above. The first is the generation of all helicity configurations and the 
second the treatment of the color summation in the case external particles with color are 
involved. 

As far as the helicity configurations are concerned this is done in a rather straight- 
forward way, resulting to an automatic evaluation of all relevant combinations. For any 
given external particle knowledge of its flavor allows the program to compute all relevant 
hehcity conflgurations. The total number is given by 2'2 3'3 where I2 is the number of 
(anti-)fermions and massless gauge bosons and Z3 the number of massive gauge bosons 
involved in the scattering process. 

The color configurations are taken into account as follows. Since only electroweak 
processes are considered in this version, the only colored particles that can appear in 
the amplitude arc quarks and antiquarks and let us have n pairs of them. Each color 
amplitude is proportional to the following color structure 

Ci = (^l,(Ti(l)(^2,o-i(2) ■ ■ ■ 5n,Ci{n) 

where CTj represents the i-th permutation of the set 1, 2, . . . , n. The code computes all 
non- vanishing color amplitudes as well as the corresponding color matrix 

and finally performs the color summation. 



3 The code HELAC 

The computation is split in two major phases. During the first phase, which we call 
initialization phase, the program selects all the relevant sub-amplitudes for the required 
process. In the simplest version the program accepts as input the following variables: 

• n the number of particles involved in the scattering 

• if 1(1 :n) the array of fiavors of the particles 
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• if lag. If set to sum over all helicity configurations is understood. If set to 1 you 
must supply also the specific helicity configuration to be computed. 

• iunitary. If set to the Feynman gauge is considered whereas if set to 1 the unitary 
gauge is used. 

• ihiggs denotes the inclusion (1) or not (O) of the Higgs particle as an intermediate 
state. 

• iwidth denotes the fixed (0) or complex (1) scheme for the introduction of the width 
of W and 

• io(l:n) is used to distinguish among initial (l) and final state particles (-1). By 
default io(l:2)=l and io(3:n)=-l. 

The first routine called after the input has been read is the routine physics/physics . f 
located in the homonymous file. In this routine all couplings of the standard electroweak 
theory are defined 0. 

Then the routine helac_init/master . f is called. In the beginning the average-over- 
helicity (avhel), the average-over-color (avcol) and the symmetry (symet) factors are 
computed. Then, by calling the routine setncc/intpar . f the number of color config- 
urations ncc is set up depending on the number of quarks and anti-quarks involved. 
For each color configuration the program constructs the skeleton of the amplitude start- 
ing at level two and proceeding up to level n — 1. All possible vertices are scanned in 
this phase, as for instance vff/panl.f which describes the coupling of a fermion anti- 
fermion to produce a vector boson, in order to select all non-zero sub-amplitudes. A 
special routine, redo/panl . f is called at the end just to check whether all the selected 
sub-amplitudes are indeed contributing to the final amplitude under consideration. The 
program ends this phase when all color configurations have been considered and the color 
matrix rmatrix(l :ncc, 1 :ncc) has been computed. 

During the second phase, which we call the computation phase, the code computes nu- 
merically the amplitude for each phase space point introduced. The main calling routine is 
helac_master/master . f . The first step is the automatic setup of the helicity configura- 
tions via the routine sethel/intpar.f . The next step is the computation of the external 
particle wave functions by the routine iniqq/panl .f . All relevant routines for this com- 
putation can be found in wavef .f. The program proceeds by computing the amplitude 
for every color configuration via nextq/panl . f . The vertex and propagator functions 
needed are included in pan2.f. Moreover the sign factor as described in Eq. (p!OD is also 

^For a recent discussion see reference (Hi. 
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computed. The routine helac_master/master.f returns the squared matrix element 
(smel2)summed and averaged over hehcity and color. Finally, a sample main/main. f 
program to run HELAC is provided. In this part of the code a call to a phase-space genera- 
tor is also included (rambo/rambo . f ), in order to generate appropriate phase space points 
for the computation. 

All floating point computations are performed in real*8 or double precision. Nev- 
ertheless the program is written in such a way that one can choose a higher accuracy if 
needed. There are two alternatives. The first one is to use real* 16 or quadruple precision. 
For this task several modules to perform quadruple precision computation with complex 
numbers are supplied in qprec.f. This is necessary because in most of the existing plat- 
forms no FORTRAN quadruple precision complex numbers are available. The most efficient 
way to use this precision with complex numbers is to exploit the virtues of FORTRAN 90 
and define the appropriate derived types and modules. The second alternative make use 
of a multi-precision library written also in FORTRAN 90. The precision is now user- 
defined and is set up by calling the routine zmset in the very beginning of the program 
(main.f). All these are automatically driven by the appropriate make files included in the 
package. 

4 Results 

In this section we will try to explain how to read the output of the program and to 
highlight several aspects of it. First of all let us follow a sample computation of the 
process 

e~e^ e~e~^'y . 

As an input we have to provide 

1. The number of particles involved in the process (5). 

2. For each particle its fiavor (2 -2 2 -2 31). In the following table we provide the 
correspondence used in HELAC: 



Ue,e , ti, (i, i/^, /i ,c, s, ... 


1,.. 


.,12 


z7g,e+,M, (i, z/^,/i+,c, s, ... 


-1,.. 


.,-12 


7,Z, W+,W- 


31,.. 


.,34 




41,.. 


.,44 



3. The parameters if lag, iunitary, ihiggs, iwidth (0 10 0). 
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The expected output is as follows 



Printed by 



Feb 03, 00 14:30 0Ut1 


Page 1/2 


Number of particles 




5 






Flavour 


of particles 




2-2 2- 


-2 31 




iflag, iunitary, ihiggs, iwidth 




10 






iunitary, ihiggs, iuni, ihi 




10 






avhel , avcol , symet 




.25 






1.0 






1.0 






for the 


1 colour conf. there are 


44 subamp 


litudes 






the number of Feynman graphs = 16 




the number of colour configurations 1 




1111 


1 




1.31323722ia50419E-04 




1111 


2 




5.70850 


J435328156E-05 




1112 


1 




.0 






1112 


2 




.0 






112 1 


1 




.0 






112 1 


2 




.0 






112 2 


1 




.0 






112 2 


2 




.0 






12 11 


1 




.0 






12 11 


2 




.0 






12 12 


1 




9.491777929411465E-06 




12 12 


2 




2.095331763482166E-05 




12 2 1 


1 





Feb 03, 00 14:30 


outi 


Page 2/2 


9 . 635699657C29954E-07 






1 2 


2 


1 


2 






1.8 


12813750492251E-06 






1 2 


2 


2 


1 






.0 












1 2 


2 


2 


2 






.0 












2 1 


1 


1 


1 






.0 












2 1 


1 


1 


2 






.0 












2 1 


1 


2 


1 






1.792073881457718E-06 






2 1 


1 


2 


2 






9.747214673799430E-07 






2 1 


2 


1 


1 






2 . 066936222894575E-05 






2 1 


2 


1 


2 






9.5924044045984e8E-0e 






2 1 


2 


2 


1 






.0 












2 1 


2 


2 


2 






.0 












2 2 


1 


1 


1 






.0 












2 2 


1 


1 


2 






.0 












2 2 


1 


2 


1 






.0 












2 2 


1 


2 


2 






.0 












2 2 


2 


1 


1 






.0 












2 2 


2 


1 


2 






.0 












2 2 


2 


2 


1 






5.708 


J69136750562E-05 






2 2 


2 


2 


2 






1.313154247256092E-04 






the 


helac amplitude is 






1 . 107657409535622E-04 







Thursday February 03, 2000 



where the first lines show the required input data; then the average-helicity, average-color 
and symmetry factors are shown; the number of contributing sub-amplitudes for each 
color configuration is given (44) as well as the total number of Feynman graphs (16). 
This marks the end of the initialization phase, and afterwards the computation of the 
squared matrix element for a given phase space point is performed. Results for each 
helicity configuration are printed and at the end the HELAC amphtude squared is given 
(1 . 107657409535622E-04) 

As far as the efficiency of the code is concerned extensive comparisons have been 
made to existing calculations. We have chosen among them two popular tools, namely 
MADGRAPH [||] and EXCALIBUR 0. Comparisons with EXCALIBUR have been restricted to 
four-fermion final states . Running under the same conditions the speed ratio was varying 
between 1 and 2. The same results have been obtained in comparisons with MADGRAPH 
taken into account that it can produce results for only up to 7 particles involved in the 
scattering process. As expected HELAC show an exponential (instead of factorial) CPU-time 
growth. There is no a priori limitation for the number of particles that HELAC can treat. 
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the only restrictions being that of memory allocation. 

In order to have a taste of a multi-precision computation we have computed the squared 
amplitude for the process 

e"(pi)e+(p2) e"(p3)e+(p4)e"(p5)e"^(P6) 

at two phase space points. Phase space point (A) is just a randomly generated one by the 
phase-space generator RAMBO. Phase space point (B) on the other hand is a very special 
one, where 

p?/GeV=100, p1+p1 = 0, p^M = 0.9, ^3 = 0, 

(P5+P6)7(P4+P5+P6)' = 0.1, ^4 = 0, (/)4 = 0, ^5 = 0, 05 = 

and nie — 0.511 x 10~^ GeV. In this case all particles are co-hnear to the beam and 
therefore important enhancements are expected. In the following table results are provided 
for these points, by using the real*8 (DP), real* 16 (QP) and 34-digit multi-precision 
version of the code. 



(A) (B) 

1.539728523150595E-008 1.256276706229023E+023 
1.53972852315058854156763002825013D-08 3.07162601093710915134136924973089D+22 
1.53972852315058854156763002825011853M-8 3.07162601093710915127950109241770808M+22 

As one can easily see the numerical stability of the DP computation is spoiled in the 
co-linear region. This is due to huge gauge cancellations occurring in this region; a way 
out of this problem will be discussed elsewhere. 



We end this presentation by summarizing the main achievements: 

• An algorithm based on Dyson-Schwinger equations has been presented that enables 
the computation of electroweak amplitudes with high efficiency. 

• Based on this algorithm a FORTRAN package HELAC has been developed, which in- 
cludes full massive computation in both the unitary and the Feynman gauge of the 
standard electroweak theory. 

• A quadruple as well as a multi-precision version of HELAC have been incorporated 
allowing numerically stable results for any phase space point and for arbitrary high 
energies. 
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Appendix 



In all calculations we are using the light-cone representation of a four-vector V^, 
defined as 



V^={V' + V,,V''-V,,V, + iVy,V,-tVy) , A=1,...A. 
Polarization state-vectors are given by 

^ f -PT PT {Px+iPy){\p\ +Pz) {Px-iPy){-\p\ + Pz 

[V2\p\' V2\p\' V2\p\pT ' V2\p\pT 

A ^ / PT -PT {Px+iPy){\p\ -Pz) {Pa, - ipy){-\p\ - Pz 



(12) 



\V2\p\'V2\p\' 



Pt 



Pt 



\P\ ^ PzPO \P\ PzPO {Px+Wy)PO iPx-iPy)PO^ 



\p\\fp^ 



(13) 



As for the Dirac matrices we are using the chiral representation 0. The wave functions 
which describe massive spinors are given by: 



«+(p) 



u-{p) 



v+{p) 



v-{p) 



( t/c \ 

aiPx + iPy)/r 
—mb/r 
V -m{px + ipy)/r J 

( m{p^ - ipy)/r \ 
—mb/r 
-a{px - ipy)/r 
V r/c j 

I -m{p^ - ipy)/r \ 
mb/r 
-a{px - iPy)/r 
\ r/c J 

( t/c \ 

a{Px + Wy)/r 
mb/r 
V m{p^ + ipy)/r J 



u+{p) 



u-{p) 



v+{p) 



v-{p) 



( mb/r \ 
m{px - ipy)/r 
—r/c 

V -a{px-ipy)/r J 

( a{px + ipy)/r \ 
—r/c 
-m{px + ipy)/r 
\ mb/r j 

( a{p^ + ipy)/r \ 
—r/c 
m{px + ipy)/r 
\ —mb/r J 

( —mb/r \ 
-m{p^ - ipy)/r 
—r/c 

V -a{p^-ipy)/r J 



(14) 



^For conventions see reference 0. 
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where: 

a = Po + \p\, b = p;, + \p\, c = 2\p\, r 
For a massless particle the spinors are 



abc 



ur{p) 



ul{p) 



I VPo + Pz \ 

iPx+iPy)/VPO+Pz 



V J 

\ 



(Px - iPy)/VPO+Pz 



Ur{p) 



V 



VPo + Pz 






-\/PO+Pz 
V -{Px - iPy)/VPO+Pz J 

( {Px + iPy)/\/P0+Pz \ 



ul{p) 



J 



V 



-^/PO+Pz 






(15) 



We now proceed to describe the vertex functions. Let us take as an example the Eq.(|D. 
By using the Dirac matrices in the chiral representation and the light-cone expression of 
four vectors, the reduced form of this equation becomes very simple, namely the four 
vector 



turns out to be 



-9R^2^J + 9Li'4'ip_l 
V -9R^l^'i + 9Li'3^2 J 



(16) 



where i = 1, . . . ,4 are the components of the spinor ip{P2) (^(A)) and 



^(1 + 75) • 



75) , ^R 

On the other hand, the spinor 

u={f+m)f{Pi)ujRiP{P2) 
used for instance in Eq.(|]), can be reduced to 

/ {-hPi + hPA)'4^i + (Ml - biPi)^2 \ 

{hP2 - &2P3)V'l + {-blP2 + hP3)i'2 

m{h2ipi - hii)2) 
\ m(-63?/'i + hii)2) I 



(17) 



u 
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where ipi, bi, pi, i = 1, . . . ,4 are the components of P, b{Pi) and ip{P2) respectively. In 
a similar way, for the standard electroweak theory in both the unitary and the Feynman 
gauge, twenty eight different vertex functions have been implemented in vertices/pan2 . f . 
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